#define A(j,i) Coeff[j*npadded+i].s0
#define SQ(a) ((a)*(a))

__kernel void QRGivensGPU_kern(
	const int npadded,
	const int n,
	__global float2* BX,
	__global float2* Coeff,
	__local float* row1,
	__local float* row2,
	__local float* B)
{
	int giX = get_global_id(0);
	int tiX = get_local_id(0);
	int ngX = get_global_size(0);
	int ntX = get_local_size(0);
	
	float c, s, tau, tau2;
	int i, j, k;
	
	if (giX < n){
		B[giX] = BX[giX].s0;	
		for(j = 0; j < n; j++) {
			row1[giX] = A(j,giX);
			for(i = n-1; i > j; i--) {
				row2[giX] = A(i,giX);
				barrier(CLK_GLOBAL_MEM_FENCE);
				if(row2[j] == 0.0f) {
					c = 1.0f;
					s = 0.0f;
				}
				else if(fabs(row2[j]) > fabs(row1[j])) {
					tau = -row1[j]/row2[j];
					s = 1.0f/sqrt(1.0f+SQ(tau));
					c = s*tau;
				}
				else {
					tau = -row2[j]/row1[j];
					c = 1.0f/sqrt(1.0f+SQ(tau));
					s = c*tau;
				}
				tau = row1[giX];
				tau2 = row2[giX];
				row1[giX] = c*tau - s*tau2;
				row2[giX] = s*tau + c*tau2;
		
		        if (giX == i || giX == j) {				
				  tau = B[j];
				  tau2 = B[i];
				}
   				barrier(CLK_GLOBAL_MEM_FENCE);
				if (giX==j) B[j] = c*tau - s*tau2;
				if (giX==i) B[i] = s*tau + c*tau2;
				
				A(i,giX) = row2[giX];
			}
			A(j,giX) = row1[giX];
		}
		barrier(CLK_GLOBAL_MEM_FENCE);
		//BX[giX].s0 = B[giX];
		for(j = n-1; j>=0; j--) {
			for(i = n-1; i > j; i--) {
				B[j] -= B[i]*A(j,i);
			}
			B[j] /= A(j,j);
		}
		//barrier(CLK_GLOBAL_MEM_FENCE);
		BX[giX].s1 = B[giX];
	}
}